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ABSTRACT 

Context. The Rosetta mission and its exquisite measurements have revived the debate on whether comets are pristine 
planetesimals or collisionally evolved objects. 

Aims. We investigate the collisional evolution experienced by the precursors of current comet nuclei during the early 
stages of the Solar System, in the context of the so-called “Nice Model”. 

Methods. We consider two environments for the collisional evolution: (1) the trans-planetary planetesimal disk, from 
the time of gas removal until the disk was dispersed by the migration of the ice giants, and (2) the dispersing disk 
during the time that the scattered disk was formed. Simulations have been performed, using different methods in the 
two cases, to find the number of destructive collisions typically experienced by a comet nucleus of 2 km radius. 

Results. In the widely accepted scenario, where the dispersal of the planetesimal disk occurred at the time of the Late 
Heavy Bombardment about 4 Gy ago, comet-sized planetesimals have a very small chance to survive against destructive 
collisions in the disk. On the extreme assumption that the disk was dispersed directly upon gas removal, there is a 
chance for a significant fraction of the planetesimals to remain intact. However, these survivors would still bear the 
marks of many non-destructive impacts. 

Conclusions. The Nice Model of Solar System evolution predicts that typical km-sized comet nuclei are predominantly 
fragments resulting from collisions experienced by larger parent bodies. An important goal for further research is to 
investigate, whether the observed properties of comet nuclei are compatible with such a collisional origin. 
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1. Introduction 

Comet nuclei are usually thought to be icy planetesimals, 
formed beyond the snow-line in the nascent Solar System. 
As such, they are naturally considered as precious targets 
of space missions - e.g ., Rosetta. This concept is supported 
by the properties of comet nuclei derived from observations. 
The low bulk densities ( e.g ., Rickman, 1989, Davidsson et 
ah, 2007), the negligible tensile strength inferred for comet 
D/Shoemaker-Levy 9 (Asphaug and Benz, 1996), and the 
low tensile strength of the surface layer required to explain 
their activity (Blum et ah, 2014) are consistent with low- 
velocity accretion, in line with the general expectation for 
planetesimals formed at large distance from the Sun. Many 
comets have exhibited evidence for a very important con¬ 
tribution by the super-volatile CO molecule to their out- 
gassing activity {e.g., Bockelee-Morvan et ah, 2004). This 
may be taken as an indication of a chemically pristine na¬ 
ture of the material that comet nuclei are made of, which 
supports the idea of a very gentle accretion process. 

However, the issue of collisional evolution in the popu¬ 
lation of icy planetesimals has also been brought to light 
in several papers. Davis and Farinclla (1997) modeled the 
collisional evolution of the Edgeworth-Kuiper Belt (EKB), 
which at the time was thought to be the source region 
for the Jupiter family comets (JFCs). They found that, 
while large EKB members are likely primordial objects, 
those with radii of a few km, like typical JFCs, are multi- 
generational fragments formed by the splitting of larger ob¬ 


jects (see also Schlichting et ah, 2013). Simultaneous with 
this first investigation, the Scattered Disk was discovered 
(Luu et ah, 1997) and it was rapidly recognized to be a 
more efficient source of JFCs than the EKB. Due to the 
longer orbital periods of its typical orbits, the Scattered 
Disk is believed to be less collisionally evolved (Rickman, 
2004), so that one could think that the observed JFCs have 
a higher chance to be primordial planetesimals than in the 
Farinella and Davis (1997) analysis. 

Another scenario was considered by Stern and 
Weissman (2001): the formation of the Oort cloud. This 
was modeled in the classical picture of gravitational ejec¬ 
tion of icy planetesimals from the growth region of the gi¬ 
ant planets (Safronov, 1977). Stern and Weissman showed 
that, during the course of this process, comet nuclei would 
be destroyed by collisions with small debris. The authors 
concluded that the storage into the Oort cloud would be 
delayed until the comet source region had been cleared of 
material so that the collisional lifetime becomes longer than 
the ejection lifetime. Naturally, in this scenario most of the 
Oort cloud comets would still bear the marks of collisional 
erosion. A similar conclusion would apply also to the origin 
of the Scattered Disk. 

Charnoz and Morbidelli (2003, 2007) - CM03/07 here¬ 
after introduced a new algorithm for evaluating the ef¬ 
fects of collisions in a population of small bodies subject to 
a complex and rapid dynamical evolution through gravita¬ 
tional perturbations, as is the case for planetesimals ejected 
from the giant planets region. This replaced the particle-in- 
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a-box models earlier used. CM03/07 showed that, for some 
appropriate initial size distributions (similar to that cur¬ 
rently observed in the trans-Neptunian region) a sufficiently 
large number of comet-size bodies would have reached the 
Oort cloud and the Scattered Disk. Although not explicitly 
discussed in these papers, they found that the vast majority 
of the Oort Cloud objects larger than 1 km in radius would 
be pristine planetesimals. However, in the Scattered Disk, 
only 2% of the final population of objects of this size would 
be primordial, the rest being collisional fragments. 

Since these papers were produced, the Nice Model for 
the early evolution of the Solar System has been introduced 
(Tsiganis et al., 2005; for the latest version, see Levison et 
ah, 2011). This changes the picture of the origin and evo¬ 
lution of comets in important ways. One central concept is 
that of the trans-planetary disk of icy planetesimals. This 
was the disk of objects formed during the infant stages of 
the Solar System beyond the original orbits of all giant 
planets, which were originally closer to the Sun. This disk 
extended out to about 30 AU and had a total mass of 20- 
50 Earth masses. It remained relatively quiescent until it 
was eventually dispersed as a consequence of a dynamical 
instability among the giant planets and of the planets’ sub¬ 
sequent migration towards their current orbits. The trans- 
planetary disk, upon its dispersal, should have given rise to 
both the Scattered Disk and the Oort cloud (Brasser and 
Morbidelli, 2013). Thus, this disk may once have been the 
repository for all the comets observed today. This would 
be compatible with the lack of evidence for any clear-cut 
differences in molecular composition (A’Hearn et al., 2012) 
or D/H ratios (Altwegg et al., 2015) between JFCs and 
comets of Oort Cloud provenance, i.e., long-period comets 
(LPCs) and Halley-type comets (HTCs). Notice that there 
is today no alternative model capable of fully explaining 
the structure of the outer Solar System without invoking a 
Nice-model-like instability of the giant planets associated 
with the dispersal of the trans-planetary disk. 

In this paper we investigate the collision rates involving 
the members of the trans-planetary disk during its whole 
evolution. First, we make the standard assumption that the 
dispersal of the disk coincided with the beginning of the 
so-called Late Heavy Bombardment (Gomes et al., 2005; 
Morbidelli et al., 2012), which means a lifetime of about 
450 Myr for the disk, prior to its dynamical dispersal. It 
is likely that the dynamical state of the disk was signifi¬ 
cantly excited. In fact, the probability that an object sur¬ 
vived the dynamical dispersal of the disk, remaining per¬ 
manently trapped into the hot EKB population (including 
mean motion resonances with Neptune) is less than 10” 3 
(see Nesvorny, 2015 for the most updated estimate); this 
means that about 1 000 Pluto-size objects should have ex¬ 
isted in the primordial disk (Stern, 1991). These bodies 
would have induced a significant excitation in the disk, 
causing a velocity dispersion of the order of 0.5 — lkm/s 
(Levison et al., 2011). If the disk stayed in this state for 
hundreds of millions of years, it is likely that the collisional 
evolution of comet-size objects has been severe. We will 
quantify this in Sect. 3.1. The conclusions will apply to 
both comets in the Scattered Disk and in the Oort Cloud, 
given that the trans-planetary disk was the progenitor of 
both these reservoirs (Brasser and Morbidelli, 2013). 

However, it is not yet certain that the dynamical disper¬ 
sal of the trans-planetary disk occurred late. The formation 
of the latest basins on the Moon (Imbrium and Orientale, 


and possibly all Nectarian basins) requires that new pro¬ 
jectiles appeared in the terrestrial planet crossing region 
several 100 My after terrestrial planet formation (Bottke et 
al., 2007). The late instability of the giant planets would 
do this in a natural way (Bottke et al., 2012; Morbidelli 
et al., 2012). However, alternative models have been pro¬ 
posed. Some (Cuk, 2012; Minton et al., 2015) invoke un¬ 
likely spectacular collision events in the inner Solar System, 
generating a large amount of debris which would have sub¬ 
sequently bombarded the Moon and the terrestrial planets. 
These models have not yet been thoroughly tested against 
all available Solar System constraints. Other models, such 
as the destabilization of a population of lunar coorbital 
objects (Cuk and Gladman, 2006) or of a fifth terrestrial 
planet that would have then partially destabilized the as¬ 
teroid belt (Chambers, 2007) did not pass such tests (Cuk 
and Gladman, 2009; Brasser and Morbidelli, 2011). The 
main argument independent of the lunar crater record in 
favor of a late dispersal of the trans-planetary disk is that 
the impact basins on Iapetus (a satellite of Saturn) have to¬ 
pographies that have relaxed by 25% or less, which argues 
that the surface layer of Iapetus was already very viscous 
at the time of basin formation. According to models of the 
thermal evolution of the satellite, this high viscosity could 
not be possible earlier than 200 My after the beginning of 
the Solar System (Robuchon et al., 2011), which implies 
that the basins of Iapetus formed late. Nevertheless, this 
constraint remains mo del-dependent. 

Thus, in the second part of the paper we consider the 
case of an early dispersal of the trans-planetary disk, oc¬ 
curring just after gas removal. In this case, the collisional 
evolution prior to the instability would be negligible (the 
disk would have survived just a few My and, presumably, 
its planetesimals would have had a very small velocity dis¬ 
persion due to gas drag); however, during the dispersal of 
the disk, the collisional evolution might have been severe 
(similar to the case studied by Stern and Weissman, 2001). 
In Sect. 3.2 we quantify the collisional evolution of comet- 
size bodies eventually stored in the Scattered Disk during 
such dispersal. Obviously, if the trans-planetary disk dis¬ 
persed late, the real collisional evolution of comets now in 
the Scattered Disk would be the sum of those studied in 
Sects. 3.1 and 3.2. 

This paper is structured as follows. We start presenting 
the logic of our reasoning, our methods, assumptions and 
choice of parameters in Sect. 2; the results are presented 
in the aforementioned Sects. 3.1 and 3.2; a discussion and 
summary of conclusions are given in Sect 4. 

2. Methods and principles 

In this work we basically follow the approach of CM03/07, 
but with some important variants detailed in this Section. 
The CM03/07 approach is very suitable to study the col¬ 
lisional evolution of a population of bodies undergoing a 
significant dynamical evolution, with orbital histories that 
can be quite different from one particle to the other. It 
basically consists of three steps. First, one does a numeri¬ 
cal simulation where the dynamical history of each particle 
is followed. Second, from the numerical simulation output, 
one computes the intrinsic collision probability and impact 
velocity of each particle with all others at each output time- 
step. Third, one assumes that each particle in the sirnula- 
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Fig. 1 . Top: the semi major axis vs. eccentricity distribution of 
the trans-planetary disk under the stirring effect of an embed¬ 
ded population of 1 000 Pluto-mass bodies, from Levison et al. 
(2011). This snapshot of the distribution is taken after 300 My 
of evolution. Bottom: the semi major axis vs. eccentricity dis¬ 
tribution of the Scattered Disk produced by the dispersal of the 
trans-planetary disk due to the giant planet instability, from 
Gomes et al. (2005). Here the snapshot of the Scattered Disk 
orbital distribution is taken 350 My after the beginning of the 
planet instability, when the Scattered Disk contains 5% of the 
original disk’s particles. 


tion is a tracer of a swarm of particles with a given initial 
size distribution. Then, using the information computed in 
the second step, one evolves the size distributions associ¬ 
ated to each tracer from one output step to another. This 
involves computing the minimal projectile size for a catas¬ 
trophic impact on targets of any given size, and the size 
distribution of the generated fragments. 

Below, we detail how we performed each of these three 
steps and in particular how we simplified the third step in 
order to reduce the parameter space we need to explore, 
although still satisfying our needs and achieving our goals. 

2.1. Numerical simulations 

In this work we use two pre-existing simulations that rep¬ 
resent well the two phases of the evolution of the trans- 
planetary disk described in the Introduction: the pre¬ 
instability phase and the dynamical dispersal phase. 


2.1.1. Pre-instability disk 

For the pre-instability phase we use one of the simulations 
of Levison et al. (2011). In those simulations, the planet 
instability occurs late and the disk is modeled in such a 
way as to mimic the self-excitation it would suffer if it 
contained 1,000 Pluto-mass bodies. As mentioned in the 
Introduction, this is a realistic number for the bodies of 
this reference mass in the original trans-planetary disk. We 
refer to Sect. 3 of Levison et al. for a technical description 
of how the simulation is done and to Fig. 2 of that paper for 
a test showing that the self-excitation process is properly 
reproduced. 

We took the state of the disk (i.e., the orbital distri¬ 
bution of the particles) 300 My after the beginning of the 
simulation. The self-stirring process increases the orbital 
excitation as \/t so that the disk is excited very rapidly, 
during the first few My, and then the evolution of the ex¬ 
citation slows down. Thus, we take the orbital distribution 
of the disk at 300 My in the simulation as representative of 
the real dynamical state of the disk during most of its pre¬ 
instability history (we will test how the results change if the 
disk’s state is taken at an earlier time in Sect. 3.1). The top 
panel of Fig. 1 shows the (a, e) distribution we consider. As 
seen, there is a clear gradient of excitation with semi ma¬ 
jor axis. This is because (i) the shorter orbital periods in 
the inner part of the disk produce more frequent encoun¬ 
ters with the massive bodies and (ii) the orbital density of 
the massive bodies is higher (in this simulation the initial 
surface density of the population of bodies in the disk is 
assumed to be proportional to 1 /r) 1 * * . 

The most up-to-date estimate for the time of the in¬ 
stability, achieved by calibrating the Nice Model on vari¬ 
ous constraints (Bottke et al., 2012; Morbidelli et al., 2012; 
Marchi et al. 2013) is ~ 4.1 Gy ago, namely 450 My after 
the disappearance of the gas from the proto-planetary disk 
(4.56 Gy ago). Given the uncertainty on this estimate, and 
to remain conservative, we assume in the following that the 
pre-instability phase of the disk lasts 400 My. 

2.1.2. The dispersal of the disk and the origin of the 
Scattered Disk 

To study the dispersal of the trans-planetary disk and the 
formation of the Scattered Disk, we use one of the simula¬ 
tions presented in Gomes et al. (2005). In that simulation 
the planets become unstable at 887 My instead of the pre¬ 
ferred date of ~ 450 My. We consider only the dynamical 
histories of the particles after the beginning of the insta¬ 
bility, ignoring the previous evolution, so that the actual 
instability date in the simulation has no influence on our 
results. In fact, the dispersal of the disk in the Nice Model 
is a very violent event, and thus the memory of what hap¬ 
pened in the pre-instability phase is quickly erased. Also, 
the overall evolution of the planetesimal population is quite 
insensitive to the exact evolution of the giant planets’ or¬ 
bits during the instability. This is evidenced by the fact 
that radically different simulations provide Scattered Disk 
populations that decay in time in very similar ways down to 

1 This assumption on the surface density profile is not arbi¬ 

trary. In fact, an accretional proto-planetary disk with a stan¬ 

dard a-prescription for its viscosity (Shakura and Sunyaev, 

1973) is expected to have a surface density profile proportional 
to r _15 A 4 (Bitsch et al., 2015). 
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~ 1% of the original disk population, as reviewed in Fig. 5 
of Brasser and Morbidelli (2013). Thus we think that the 
simulation that we consider is sufficient for our purposes. 

The simulation covers a timespan of 350 My after the in¬ 
stability, identifies each particle individually and produces 
a Scattered Disk made of 5% of the original particles at 
this date, whose (a, e) distribution is depicted in the bot¬ 
tom panel of Fig. 1. 


2.2. Collision probabilities and velocities 

The simulations we consider record the orbital elements of 
all particles at regular output intervals dt ou t- Given any pair 
of particles, we compute their intrinsic collision probability 
and impact velocity using the Opik-like algorithm described 
in Wetherill (1967), implemented in a fortran code by P. 
Farinella and D. Davis and kindly provided to us. The algo¬ 
rithm considers the orbital elements a, e, i (semi major axis, 
eccentricity and inclination) of each of the two particles and 
assumes that the orbital angles w,D,M (argument of per¬ 
ihelion, longitude of node, mean anomaly) precess linearly 
with time (so that their values are random on a sufficiently 
long time interval), without inducing changes on (a,e,i). 
Clearly, these are approximations, but they are basically 
correct until the particles reach very large inclinations and 
undergo large amplitude Kozai cycles (Vokrouhlicky et al., 
2012; Pokorny and Vokrouhlicky, 2013), which is not the 
case for the pre-instability disk nor for most of the particles 
in the Scattered Disk (Kozai cycles for trans-planetary or¬ 
bits are pronounced only in mean motion resonances or for 
planet-crossing orbits; Thomas and Morbidelli, 1996). The 
use of a collisional probability algorithm like Wetherill’s 
on the output of a numerical integration is standard prac¬ 
tice and leads to quite accurate results ( e.g ., Levison et al., 
2000; Rickman et al., 2014). 

The code returns the intrinsic collision probability Pi, 
which is the probability that a point-like projectile hits a 
R = 1 km target in an year. Thus, the probability that two 
objects of radii Ri and R 2 (expressed in km) collide over 
a time interval 5t is therefore P co u = Pi(Ri + R 2 ) 2 5t. The 
impact velocity v co u is the mean of the relative velocities 
between the two orbits over all collision configurations. It 
corresponds to the velocity of approach, prior to any ac¬ 
celeration due to the mutual attraction between the two 
bodies. The latter is negligible for planetesimals. 

For the pre-instability disk, it is not necessary to con¬ 
sider the collision probability of each of the simulated parti¬ 
cles. Averaged values are enough. However, given the radial 
excitation gradient shown in the top panel of Fig. 1, we di¬ 
vide the disk in three zones: zone I with 15 < a < 20 AU, 
zone II with 20 < a < 25 AU and zone III with 25 < a < 
30 AU. Then, denoting by k the particles in one zone and 
to those in the other zone (possibly the same zones), we 
compute the mean intrinsic collision probability as: 

K M 

^ = ’ (1) 
k =1 m— 1 

where K and M are the total numbers of particles in the 
considered disk zones and Pi(k, to) is the intrinsic probabil¬ 
ity between particles k and to. Similarly, the mean impact 


velocity (weighted by collision probability) is: 

1 KM 

Vc°u = KM p. (k, to) , (2) 

* fc=l m= 1 

where v co u(k,m ) is the collision velocity between particles 
k and to. 

For the simulation of the disk dispersal, instead, we com¬ 
pute the collision probability individually for each particle 
that will be a dynamical survivor in the Scattered Disk 
at the end of the simulation, against all other particles. 
Denoting by j a Scattered Disk particle and by l any other 
particle, the mean collision probability of particle j at time 
t is therefore: 


1 Ht) 

PiU,t) = J^Y, p iU, 1 ) > (3) 


where L(t) is the total number of particles surviving in the 
integration at time t. For the velocity, we have: 


L(t) 


Vcoll{j,t) = 


L{t)Pi(j,t) ^ 


'Y^P i {j,l)v co u{j,l) • (4) 


2.3. Size distributions and disruption probabilities 

In the method introduced in CM03/07, one defines an ini¬ 
tial size distribution for the swarm of planetesimals repre¬ 
sented by each simulation particle. At each step, using the 
pre-computed collision probabilities among the simulation 
particles, one computes the number of collisions occurring 
between planetesimals of any given size. From the impact 
velocities and masses of projectile and target, one computes 
the consequences of the collisions (cratering event, catas¬ 
trophic break-up) and the size distribution of the generated 
fragments. In this way, the evolution of the size distribu¬ 
tions associated to each simulation particle is computed. In 
the end, one evaluates which fraction of the planetesimals 
of a given size are survivors of the original population or 
collisional fragments of larger planetesimals. 

This approach is correct, but it is overshooting for our 
goal in this paper, which is just to assess whether a comet- 
size object is likely to have avoided catastrophic collisions. 
Moreover, it requires exploring a variety of initial size dis¬ 
tributions, demanding a quite tedious exploration of the 
parameter space defining them. 

Thus, we modify and simplify the approach as described 
below. We use the approach typical of a mathematical 
demonstration ad absurdum (by reduction to the absurd). 
That is, we start by assuming that the planetesimals down 
to comet-size objects are not significantly affected by col¬ 
lisions. This means that the planetesimal size distribution 
does not evolve with time, and that the initial distribu¬ 
tion has to be the same as the current distribution in the 
Scattered Disk, just scaled up by the inverse of the implan¬ 
tation efficiency (the fraction of the disk population surviv¬ 
ing in the end in the Scattered Disk). We detail below what 
this size distribution is. 

Then, using this distribution and the pre-computed col¬ 
lision probabilities and velocities, we evaluate the mini¬ 
mum size of a projectile capable of disrupting a comet-size 
body and thus the number of catastrophic collisions n co u 
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that each comet-size body should suffer (we detail below 
how we do this evaluation). The probability that a comet- 
size body has escaped all catastrophic collisions is then 
Pintact = exp(-n coJ z). If Pintact is close to unity, then our 
assumption of a negligible collisional role is verified. But if 
Pintact is small, then we reach the absurd situation that, by 
assuming that the planetesimal population was not affected 
by collisions, we conclude that most planetesimals should 
have been destroyed! This means that the assumption was 
wrong, and hence, that the planetesimal size distribution 
was significantly affected by collisions. 

With this approach we cannot compute the actual prob¬ 
ability that a comet-like body has escaped all catastrophic 
collisions, but we know that it has to be smaller than 
Pintact- In fact, any initial size distribution evolving by col¬ 
lisions towards the current Scattered Disk distribution must 
have had originally more bodies than we assumed (because 
of collisional comminution) and therefore, the probability 
that a given body was catastrophically disrupted must be 
higher than we computed (because of a larger initial num¬ 
ber of projectiles). If the value of Pintact that we computed 
is already small, this is enough for our purposes. 

In this paper, as comet-size bodies we consider ob¬ 
jects with radius R = 2 km, appropriate for comet 
67P/Churyumov-Gerasimenko, the target of the Rosetta 
mission (see Rickman et ah, 2015). 

2.3.1. The Scattered Disk size distribution and the minimal 
number of comet-size objects in the original 
trans-planetary disk 

The most recent estimate of the Scattered Disk population 
has been presented in Brasser and Morbidelli (2013). In 
that work, as in Duncan and Levison (1997), the number of 
comet-size bodies in the Scattered Disk is evaluated from 
(a) the number of known Jupiter family comets in some 
given range of orbits and magnitudes for which the JFC 
sample is assumed complete and (b) the numerical relation¬ 
ship between the Scattered Disk population and the Jupiter 
family population that the former sustains, obtained from 
numerical simulations. With respect to previous estimates 
(e.g., Duncan and Levison, 1997), the estimate in Brasser 
and Morbidelli is improved in two respects: it uses the most 
recent conversion from total magnitude to nuclear size from 
Fernandez and Sosa (2012) and it is based on new sim¬ 
ulations deriving Jupiter-family comets from a Scattered 
Disk that is excited in inclination (the original Duncan and 
Levison work assumed that inclinations in the Scattered 
Disk are of a few degrees only, which has then been re¬ 
futed by observations). Brasser and Morbidelli concluded 
that there are 2 x 10 9 bodies in the Scattered Disk to¬ 
day larger than 2.3 km in diameter. Given a Scattered Disk 
implantation efficiency of 1%, this means that the original 
trans-planetary disk, if not affected by collisional comminu¬ 
tion, should have contained 2 x 10 11 of these bodies. 

We believe that this is a lower bound for the original 
disk population for three reasons. First, in a subsequent 
work accounting for an improved fading law (probability 
that a comet does not survive more than n perihelion pas¬ 
sages), Brasser and Wang (2014) raised the estimate of the 
Scattered Disk population to 6 x 10 9 bodies larger than 
2.3 km in diameter. Second, serendipitous stellar occulta- 
tion observations by the HST guiding sensors (Schlichting 
et ah, 2009) and by the Corot survey (Liu et ah, 2015) 


suggest that the average sky density of bodies larger than 
250m in radius over a ±5° ecliptic band is 2 x 10 7 /deg 2 . 
Thus, there are at least 7 x 10 10 bodies of this size in the 
trans-Neptunian region; with a cumulative size distribu¬ 
tion proportional to R~ 2 this implies 3.5 x 10 9 objects with 
D > 2.3 km. It is unclear which population (cold EKB, 
hot EKB, Scattered Disk) the detected objects belong to. 
But, given that the Scattered Disk outnumbers the others 
(compare Trujillo et al., 2000 with Fraser et al., 2014 for 
the observational point of view and Brasser and Morbidelli 
2013 with Nesvorny, 2015 for the modeling point of view), 
the number above can be considered to be an estimate - 
if not a lower bound - of the Scattered Disk population. 
Third, repeating the same exercise for the Oort cloud pop¬ 
ulation, Brasser and Morbidelli (2013) estimated that the 
primordial trans-planetary disk should have contained 10 12 
objects with D > 2.3 km; thus probably the reality lies in 
between 2 x 10 11 and 10 12 . 

Thus, to remain conservative (i.e., underestimate the 
total number of collisions), in this work we assume that 
the trans-planetary disk contained 2 x 10 11 objects with 
D > 2.3 km. As for the size distribution, we turn again 
to comet observations. Estimates of the JFC size distri¬ 
bution range significantly from authors to authors, from 
quite steep (exponent of the differential distribution close 
to —3.5 - Fernandez et al., 1999; Tancredi et al., 2006 - or 
even steeper - Belton, 2015) to shallow (differential slope 
of —2.6; Lowry and Weissman, 2003). Consequently, in this 
work we assume as nominal differential slope the value —3 
(Meech et al., 2004; Lamy et al., 2004; Snodgrass et al., 
2011), but we also study the dependence of the results 
on exponents for the differential distribution ranging from 
-2.5 to -3.5. 

A shallow size distribution is preferred according the the 
most recent planetesimal formation models (Johansen et 
al., 2015: q = —2.8). TNO surveys (e.g. Fraser et al., 2014) 
also suggest that the size distribution of objects smaller 
than 50 km in radius is shallow (q between —3.1 and —2.5, 
although it may steepen up for not yet detectable comet- 
size bodies). 

For reference, a disk with a size distribution similar to 
that of the hot EKB (Fraser et al., 2014), namely with a dif¬ 
ferential slope of —3 for R < 50 km and —5 for R > 50 km, 
a total number of 2 x 10 11 objects with R > 1.15 km and 
a density of 1 g/ cm 3 would have a total mass of 35 Earth 
masses, in good agreement with that required by the Nice 
Model (Gomes et al, 2005; Morbidelli et al., 2007; Nesvorny 
and Morbidelli, 2012). 

2.3.2. Minimum size of catastrophic projectiles 

The kinetic energy of an impact that can catastrophically 
destroy an object is 

E disr upt = ^npR 3 Q*(R) , (5) 

where p is the bulk density of the target of radius R; Q*{R) 
is the specific energy for disruption and is size dependent. 
There are several Q* ( R ) laws proposed in the literature for 
various materials. Benz and Asphaug (1999) propose two 
such laws for bodies made of ” strong ice”, hit respectively 
at 1 km/s and 3km/s. As we will see in Sect. 3, the former 
velocity is well adapted to the pre-instability disk while the 
second is suitable for the disk dispersal phase. For a R = 
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2 km body the two values of Q*(R) are actually comparable. 
Leinhardt and Stewart (2009) produced a Q* law for bodies 
made of ’’weak ice”, hit at lkm/s. Their Q*(R) function 
follows the general trend of those in Benz and Asphaug 
(1999) but the value for R = 2 km is about an order of 
magnitude smaller (see their Fig. 11). The scaling of Q* 
with velocity given in Eq. (2) of Stewart and Leinhardt 
(2009) gives a value 2.4 times larger for a velocity of 3 km/ s, 
which is still 4 times smaller than that reported in Benz and 
Asphaug (1999) for the same speed. Even if the Leinhardt 
and Stewart value may be more appropriate for pristine, 
low-density planetesimals, we use the Benz and Asphaug 
values in order to be conservative once again. This likely 
overestimates the minimal size of a projectile capable of 
disrupting the target and thus underestimates the number 
of catastrophic collisions. 

In fact, once Edisrupt is known, the minimal size of a 
catastrophic projectile r p is given by the equation 

2 n P r p2 V coll — Edisrupt j ( 6 ) 

so that, the higher is Edisrupt the larger is r p . Notice that, 
if one assumes that the bulk density of projectile and target 
is the same, p simplifies from the right hand and left hand 
sides of (6) and the result is independent of p. 


Table 1 . Table of results for the pre-instability disk. The 
first row reports the disk zone where the target is located 
and the first column reports the disk zone where the pro¬ 
jectile is located. The disk zones are: (I) a < 20 AU, (II) 
20 < a < 25 AU, (III) a > 25 AU. Then, each box re¬ 
ports on the top the mean intrinsic collision probability Pi 
(number of collisions per year per projectile on a target of 
R = 1km), in the middle the mean collision velocity v co u 
(in krn/s) and at the bottom the minimum size of a catas¬ 
trophic projectile r p (in km). 


\ target 
projectile \ 

i 

II 

III 

i 

1.85 x 10~ 20 
0.78 

0.23 

3.75 x 10" 21 
0.74 

0.24 

1.00 x 10 -24 
0.95 

0.20 

ii 

3.75 x 10“ 21 
0.74 

0.24 

8.95 x 10" 21 
0.44 

0.33 

7.95 x 10“ 22 
0.38 

0.37 

in 

1.00 x 10“ 24 
0.95 

0.20 

7.95 x 10“ 22 
0.38 

0.37 

7.32 x 10 -21 
0.24 

0.51 


2.3.3. Total number of catastrophic events 

Once the minimum size of a catastrophic projectile is 
known, the total number of catastrophic impacts for a tar¬ 
get of radius Rt is computed as: 

r Rmax 

Neon = (Pi5t) / ( R t + R p ) 2 N(R p )dR p , (7) 

* r p 

where Pi is the considered intrinsic probability (averaged 
over the ensemble of potential projectiles, as explained in 
Sect. 2.2), St is the considered time-span, N(R p )dR p is the 
differential size distribution, r p is the minimum size for 
a catastrophic projectile and R m ax is the maximum size 
for which the considered size distribution is valid. Given 
that the size distribution of the trans-Neptunian popula¬ 
tions turns from steep (at the large size end) to shallow (at 
the small size end) at a size of approximately R ~ 50 km 
(Bernstein et al., 2004; Fuentes et al., 2009; Fraser et al., 
2014), we assume Rmax = 50km. We neglect the relatively 
small contribution by projectiles of even larger sizes. 

Eq. (7) is often approximated by 

r R ma* 

Neon = PiStRT / N(R p )dR p = PiStR^N(> r p ) , 

J r p 

( 8 ) 

where N{> r p ) is the cumulative number of bodies larger 
than r p . This approximation is good for steep size distri¬ 
butions, or in the limit r p —> 0. However, for shallow size 
distributions like the one we consider here and r p not much 
smaller than Rt (as is the case for low velocity collisions), 
the approximation is not very precise. Thus, we solve the 
integral (7) exactly. That is, denoting by q the exponent of 
the differential size distribution, the primitive of the inte¬ 


grand in (7) is: 

q = -2.5 : -2 (R^ + 6 R T R P - 3i?2)/(3i?p /2 ) 

q = -3 : ~R^/(2R^) - 2 R T /R P + log (Rp) 

q = -3.5 : -2(3 R%, + 10 R T R P + 15Rl)/(15R 5 p /2 ) . 

( 9 ) 

3. Results 

We report here the results obtained for the pre-instability 
disk and the disk dispersal phase, obtained applying the 
methods described in the previous Section. 

3.1. Disruptive collisions in the pre-instability disk 

We show in Table 1 the results concerning the intrinsic 
collision probability Pi, the collision velocity v co u and the 
minimum size of a catastrophic projectile r p for a target 
with R = 2 km, considering all possible combinations of 
disk zones for projectile and target. 

As seen, the intrinsic collision probability is higher if 
both target and projectile are in the inner part of the disk 
than in the outer part. The collision velocity is also higher. 
The mean intrinsic collision probability is lower if target 
and projectile belong to different disk zones, because not 
all particle orbits from the two zones intersect. 

Table 2 reports the number of catastrophic collisions 
expected for a 2 km target, for the three considered values 
of the exponent q of the differential size distribution. The 
calculation has been done assuming St = 400 My (the ex¬ 
pected lifetime of the pre-instability phase), and applying 
(7) and (9) to the numbers reported in Table 1 for projec¬ 
tiles in each disk zone. Because the surface density of the 
disk is assumed to be proportional to 1 /r in the Levison et 
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Table 2. Number of disruptive collisions expected for a 
R = 2 km target located in each disk zone, as a function 
of the exponent q of the differential size distribution. The 
first row reports the target’s zone. The first column gives 
the value of q. Each box reports the number of catastrophic 
collisions expected over 400 My. In parentheses we report 
the same quantity estimated by using the dynamical state 
of the disk after 100 My of evolution, instead of that shown 
in the top panel of Fig. 1 (300My). The number of catas¬ 
trophic collisions is smaller, but it is nevertheless much 
larger than unity in all cases. 


^ytargetzone 

I 

II 

III 

-2.5 

58.0 (51.2) 

28.7 (20.7) 

12.3 (9.6) 

-3.0 

94.5 (75.0) 

39.7 (23.7) 

12.1 (7.9) 

-3.5 

190.6 (137.7) 

70.2 (35.3) 

15.4 (8.2) 


Table 3. The same as Table 2 but now reporting the radius 
of a body (in km) for which the probability of catastrophic 
impact is 50%, for each disk zone and assumed slope of the 
projectile size distribution. In parentheses we report the 
same quantity estimated by using the dynamical state of 
the disk after 100 My of evolution, instead of that shown in 
the top panel of Fig. 1 (300 My). When the result exceeds 
50 km (the value at which the size distribution is no longer 
described by a power-law with exponent q), we report 50 
for simplicity. 


^^targetzone 

I 

II 

III 

-2.5 

50 (50) 

50 (50) 

50 (43) 

-3.0 

50 (50) 

50 (50) 

48 (37) 

-3.5 

50 (50) 

50 (44) 

29 (20) 


al. (2011) simulation that we use, an equal number of pro¬ 
jectiles of a given size is assumed to exist initially in each 
of the disk zones. 

We notice that the total number of collisions for comets 
in zone III of the disk has a very weak dependence on q , 
because the size of the minimum catastrophic projectile is 
quite big (0.5km in radius). The opposite is true for comet- 
size targets in zone I of the disk. Clearly, in all cases the to¬ 
tal number of collisions is larger than 1. This means that the 
probability that a 2 km body escapes from all catastrophic 
collisions is small. From the numbers in Table 2 this prob¬ 
ability is always smaller than exp(—12) = 6 x 10“ 6 ; using 
the numbers in parentheses, we get exp(—7.9) = 4 x 10 -4 . 

Table 3 reports the radius of the comets for which the 
probability of a catastrophic impact over the age of the disk 
is 50%, again for each disk zone and assumed value of q. 
This size is extremely large, in most cases exceeding 50 km. 
This is beacuse the number of catastrophic impacts depends 
weakly on the size of the target. In fact, if one assumes 
for simplicity that Q* is independent of size, the size of 
the minimal catastrophic projectile scales linearly with the 
target size Rt\ then, using (8) one finds that the number 
of catastrophic impacts scales as R^~ q+1 which, for q = 
—3 eliminates the dependence on Rt . For the compilation 
of Table 3 we nevertheless used the dependence of Q* on 


radius given in Benz and Asphaug (1999) and the non- 
approximated formulae (9). 

Note that this result is valid for both Scattered Disk 
comets (JFCs) and Oort cloud comets (LPCs/HTCs), be¬ 
cause both reservoirs form from the same disk (Brasser and 
Morbidelli, 2013). Thus we conclude that if the giant planet 
instability occurred as late as 4.1 Gy ago, the possibility 
that a 2 km comet is a pristine planetesimal, rather than a 
collisional fragment, is very slim. 

3.2. Disruptive collisions during disk dispersal 

As we said in Sect. 2.2, for the disk dispersal phase we 
consider individually each particle ending in the Scattered 
Disk, because their orbital histories can be very diverse. 
However, because the assumption in our ad absurdum ap¬ 
proach is that the disk is not collisionally active and its size 
distribution does not evolve, for each target particle j we 
can average over time the value of Pi(j,t) given in (3) as: 

= ( 10 ) 

and apply the result over a time interval St = 350 My, which 
is the integration timespan. 

However, this 350 My simulation timespan covers only 
a small fraction of the lifetime of the Solar System and in 
principle there may still be a significant collisional evolu¬ 
tion in the Scattered Disk over the remaining ~ 4 Gy of 
Solar System history. To have an estimate of the collision 
probability over this remaining time, we proceed as follows. 
We assume that the orbital distribution in the Scattered 
Disk does not evolve with time, remaining equal to that 
shown in the bottom panel of Fig. 1, but we assume that the 
Scattered Disk population decays with time. The Scattered 
Disk at the end of the simulation accounts for 5% of the ini¬ 
tial disk particles and over the rest of the Solar System life¬ 
time it would decay to about 1% (Brasser and Morbidelli, 
2013). So, we assume that number of Scattered Disk par¬ 
ticles decays as exp (—t/r) with r such that after 4 Gy the 
population is reduced by a factor of 5. The integral of the 
exponential function over 4 Gy with such a value of r is 0.5. 
Thus, we take the last computed value of P(j,t) for each 
Scattered Disk particle (t = 350 My) and we multiply it 
by St = 4 Gy and then divide by 2. We find that the inte¬ 
grated collision probability over the last 4 Gy is a fraction 
(typically from 10% to 80%) of that integrated over the 
first 350 My. This is because at the beginning the disk is 
much more populated and, moreover, the most collisional 
active phase is when the disk is just stirred up by the plan¬ 
ets’ action (Stern and Weissman, 2001). Instead, once the 
Scattered Disk is formed, the collision probability per unit 
time per particle is strongly reduced, due to the large or¬ 
bital space that the Scattered Disk fills and the long orbital 
periods. 

Given the typical collision velocities of 2 — 4 km/ s, the 
typical size of the smallest catastrophic projectile for a 
R = 2 km target is ~ 100 m. One may wonder whether 
bodies this small existed in the disk. The occultation obser¬ 
vations mentioned above (Schlichting et ah, 2009; Liu et ah, 
2015) show that bodies of this size exist today. Following 
our ad absurdum approach, we then need to assume that 
they existed in the original disk, because if the disk did not 
evolve collisionally, no objects could be generated. 
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Fig. 2. The number of expected catastrophic collision for each 
particle surviving in the Scattered Disk at the end of the disk 
dispersal simulation. The symbols depict different values for the 
exponent of the differential size distribution q, as labeled in the 
plot. 

Fig. 2 shows for each Scattered Disk particle the ex¬ 
pected number of catastrophic collisions that a 2 km target 
should have suffered; the colors refer to different values of q. 
The number of catastrophic collisions changes considerably 
from particle to particle, because the dynamical histories of 
Scattered Disk objects can be very diverse. We see that, if 
q = —3.5 or steeper, clearly each comet-size object should 
have suffered at least two catastrophic collisions, with an 
average of 4.7 collisions. The radius of comets for which 
the average number of catastrophic impacts would be 0.5 
is approximately 10 km. This excludes the possibility, pro¬ 
posed by Belton (2015), that the size distribution of the 
Scattered Disk is steeper than q = —3.5 below this radius. 
In fact, the disk would be collisionally evolved and there¬ 
fore it would have acquired a collisional equilibrium size 
distribution, which implies q = —3.5 (Dohnanyi, 1969) or 
shallower (|q| < 3.5; O’Brien and Greenberg, 2003). 

If the distribution is very shallow (q = —2.5), Fig. 2 
shows that about half of the comets should have had no 
catastrophic collisions, the average number of catastrophic 
collisions per object being 0.5. The nominal case q = —3.0 
is borderline. Most comets should have had at least one 
catastrophic collision (the average being 1.5 collisions per 
comet) but some, with favorable orbital histories, would 
have had no collision at all. Please notice that, had we used 
the 4 times smaller value of Q* from Leinhardt and Stewart 
(see Sect. 2.3.2), the number of catastrophic impacts would 
have increased by a factor of ~ 3 for q = —3.5, a factor 
~ 2.5 for q = —3 and a factor of 2 for q = —2.5. 

Fig. 3 shows the same results but using a different repre¬ 
sentation. The number of collisions N co u is converted into 
a probability to avoid all collisions P(0) = exp (—N co u). 
Then, the normalized cumulative distribution of the P(0) 
values is plotted. The thick curves are for the nominal Q* 
value from Benz and Asphaug (1999) and the thin curves 
for a Q* value 4 times smaller. 

Note that, if the disk dispersal occurred late, the num¬ 
ber of catastrophic collisions found in this Section should be 
added to those reported in Table 2. So, all original comets 
should have been destroyed. If instead the dispersal of the 
disk occurred soon after the disappearance of gas from the 
disk, the results illustrated in Figs. 2 and 3 apply alone. In 
this case, if the original and current size distributions in the 



Fig. 3. Fraction of particles ending in the Scattered Disk with 
a probability to escape all catastrophic collisions P(0) smaller 
than that indicated on the horizontal axis. This is an alterna¬ 
tive represenation of the results already shown in Fig. 2. The 
different line styles refer to different exponents for the differen¬ 
tial size distribution q, as labeled on the plot. The thick curves 
correspond to the value of Q* given in Benz and Asphaug (1999) 
and the thin curves to a Q* value 4 times smaller, as in Leinhardt 
and Stewart (2009) with the velocity scaling provided in Stewart 
and Leinhardt (2009). 

comet-size range are shallow and a quite large value of Q* 
applies, there is the possibility that a fraction of the comets 
are pristine objects which escaped catastrophic collisions. 

4. Conclusions 

In this work we have estimated the total number of catas¬ 
trophic collisions that a typical Jupiter family comet (here 
assumed to have radius R = 2 km) should have had over 
its dynamical lifetime, first in the trans-planetary disk and 
then in the Scattered Disk. 

We have shown that, if the trans-planetary disk be¬ 
yond the original orbit of Neptune has been dispersed by 
a late dynamical instability of the giant planets occurring 
~ 4.1 Gy ago, comet-size objects should have suffered nu¬ 
merous catastrophic collisions in the pre-instability phase. 
Thus, not only JFCs, but also Oort cloud comets should 
be fragments of originally larger bodies. Because the late 
instability of the giant planet system is, at the current level 
of understanding, the best explanation for the trigger of the 
Late Heavy Bombardment of the Solar System, the forma¬ 
tion of late lunar basins on the Moon (Bottke et ah, 2012; 
Morbidelli et ah, 2012) and the impact age record on me¬ 
teorites (Marchi et al., 2013), we believe that this is the 
conclusion of our work. 

However, in the hypothetical case that the dispersal of 
the disk occurred early, the collisional evolution of comet- 
size bodies ending in the Scattered Disk would have been 
less severe. If the size distribution of comet-size objects 
in today’s Scattered Disk and in the primordial trans- 
planetary disk was shallow (differential index |qj < 3), it 
is possible in principle that a significant fraction of comet- 
size objects escaped all catastrophic collisions. 

The reader should keep in mind, though, that through¬ 
out our study we have taken the most conservative assump¬ 
tions, so that the number of catastrophic collisions that 
we computed should be considered as a lower estimate. 
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In fact, we have considered an initial size distribution in 
the disk that contains the minimum possible number of 
comet-size objects (Sect. 2.3.1). Also, we have assumed the 
specific energy for catastrophic disruption given in Benz 
and Asphaug (1999), which probably overestimates that ap¬ 
propriate for weak icy aggregates (Leinhardt and Stewart, 
2009; Sect. 2.3.2): the adoption of the 4 times smaller spe¬ 
cific energy for disruption of Leinhardt and Stewart (2009) 
would have multiplied the number of catastrophic impacts 
shown in Fig. 2 by a factor 2.5 for q = —3 and a factor of 
2 if q = —2.5, while producing the thin curves in Fig. 3. 
Moreover, the ad absurdum approach that we followed by 
its essence provides just a minimal estimate of the num¬ 
ber of collisions (Sect. 2.3). Finally, one should take into 
account that, for each catastrophic collision, the number 
of quasi-catastrophic collisions would be much higher (be¬ 
ing caused by smaller projectiles, which are more numer¬ 
ous). Thus, even if a comet had not suffered, by chance, 
any catastrophic collision, its morphology would have been 
sculpted but numerous large sub-catastrophic impacts. 

Therefore, we conclude that typical JFCs of the size 
of 67P/Churyumov-Gerasimenko most likely are not in¬ 
tact planetesimals but they are either fragments of origi¬ 
nally bigger bodies (the most likely case), or are planetesi¬ 
mals strongly sculpted by barely sub-catastrophic impacts 2 . 
In the latter case, this sculpting may disappear after the 
comets have been severely eroded by sublimation or split¬ 
ting. However, some JFCs are statistically likely to bear 
these scars. 

There are some properties of comets that appear to 
be in contradiction to a collisional origin and would 
rather suggest that comets are primordial survivors. In the 
Introduction we mentioned the low bulk densities and neg¬ 
ligible tensile strengths along with the large abundance of 
super-volatiles like CO. These properties are certainly com¬ 
patible with an origin of comets as primordial survivors 
of the icy planetesimal population. Hence, our results lead 
to the obvious question of whether collisions are able to 
conserve these primordial properties in the fragments they 
produce. This question is currently unanswered and merits 
careful consideration. 

If a detailed modeling of the collisions between icy plan¬ 
etesimals would show that the primordial-like features of 
comets are not preserved in the fragments, one may suspect 
that our current vision of outer Solar System evolution is 
not appropriate. For instance, there might have been no de¬ 
lay in the dynamical instability, or the disk remained less 
self-excited due to a smaller number of large bodies than 
we envision, or there was a drastic cut-off in the size distri¬ 
bution affecting sub-km objects, thus limiting the number 
of catastrophic projectiles. 
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